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SUMfiARY ‘ ‘ 

A multi-level grid method has been studied as a possible means of 
accelerating convergence in relaxation calculations for transonic flo»vs. The 
method employs a hierarchy of grids, ranging from very coarse (e.g. 8x2 
mesh cells) to fine (e.g. 128 x 32); the coarser grids are used to diminish 
the magnitude of the smooth part of the residuals, hopefully with far less 

total work than would be required with, say, optimal SLOR iterations on the 
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finest grid. The method was applied to the solution of the transonic small- 
disturbance equation for the velocity potential in conservation fonn. 

Nonlifting transonic flow past a parabolic-arc airfoil is the example studied, 
with meshes of both constant and variable step si7.e, 

INTRODUCTION 

The multi-level grid method, for accelerating convergence in relaxation 
calculations, has been shown to be very efficient for solving elliptic problems 
with Oirichlet boundary conditions. For background and historical material, 
see references 1 to 4. In reference 5, Brandt gives an extensive discussion 
and analysis of the method, together with several different procedures for 
applying the method, The idea of the method is based on the fact that in many 
typical elliptic boundary-value problems, the err.or is composed of a discrete 
spectrum of wave lengths, which range froi« the width of the region down to the 
width of a mesh cell. The short wave-length components of the error are 
usually diminished quite rapidly in a relaxation calculation, while the long 
wave-length components diminish very slowly, After only a few iterations the 
residual will be smooth, since the short wave-length error components have 
been eliminated; and thus the residual can be represented accurately on a 
coarser mesh. An equation called the "residual" equation is then solved on 
the coarser mesh, and the resulting correction is added to the last approxi- 
mation on the fine mesh, yielding a significant improvement with very little 
work. 

Since relaxation methods are currently the most attractive for obtaining 
numerical solutions to transonic aerodynamics problens, the question arises 
as to whether a multi-level, or multi-grid (MG) method can be used in a mixed 


E 


flow with shock waves. In this paper we report some early results using the 
MG method to solve a simple transonic problem: we consider the tran- 
sonic small-disturbance equation for the velocity potential, for nonlifting 
flow past a parabolic-arc airfoil. 

PROBLEM DESCRIPTION 


The transonic small -disturbance equation for the velocity potential can 
be written in conservation form as: 

Px + qy = 0 (1) 


where 


■[ 


K - (x+ 1) m£ * 


'00 




( 2 ) 


q " ^y (3) 

K = (1-M2)/t2/3 ‘ (4) 

Equation (1) is to be solved subject to the boundary conditions that the 
disturbance potential, (J), vanishes at infinity and the flow is tangent to 


the airfoil surface, in the interval 


^ 1/2 ; i .e. 


at y = 0, <|)y = F'(x) for 
= 0 for , 


1/2 


> 1/2 

ace) thickness distribution function. 


(5) 


where F(x) is the (upper suri 
t is the usual thickness ratio, and Y, Mm, and K are the ratio of specific 
heats, free-stream Mach number, and transonic similarity parameter, respec- 
tively. The form of equations (1) to (5) is a correctly-scaled transonic 
similarity form, in that all quantities are of order 1. Physical quantities, 
denoted by a "hat" symbol are related to the scaled quantities as follows: 
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( 6 ) 


X “ CX 

J s 

t (x) ^ 2 ctF(x) 

where c 1s the airfoil chord length and t is the total thickness distribution 
of the syninietric airfoil. 

Equation (1) is of hyperbolic or elliptic type depending on whether 

U = K-(y + l)M^<})x (7) 

is negative or positive, respectively, 

Finite-Difference Equations 

Munnan's conservative difference scheme (ref. 6) can be conveniently 
presented in teniis of Jameson's "switching function" (ref. 7) as follows: 


(1 - p..) P.. + p . . P. . . + Q. = 0 


where 


Pij = u;. 


. ♦jil ,j ~^*ij 


Ax 


U,.J = K -(y+i)m£ 


« ♦i.-n-l 

^i I = 2 


( 8 ) 

(9) 

( 10 ) 


01) 


( 12 ) 
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and where 





= 1 


if U... > 0 

. , 0 


02 ) 


It should be ooted here that, in the interest of simplicity, we have 
presented only the constant-step-size (unstretched grid) fonii of the dif- 
ference equations. In the case of a stretched grid, the conservative 
difference equations cannot be factored into the nice form given above, but 
this presents no real difficulty. The actual computer program is written 
for a stretched grid, with the identity transformation (constant step size) 
included as a special case. 

Vertical Line Relaxation 

A vertical line relaxation scheme for solving ’equation (8) by iteration 
can be written as: 






(13) 


Where T, j = ♦ ij - ♦ 1j O'*) 

if,’*' denotes a "new" value of <}), obtained during the latest iteration sweep, 
while <p is the value from the previous sweep. R. ., which is the left-hand 

• J 

side of equation (8), is evaluated with "old" values of (f. ., as are the 

* J 

iteration coefficients A through E, which are given in the appendix. 
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Multi-Grid Approach 

Let us introduce a sequence of grids G^, Gg, where for simplicity, 

" ^^k+1, >'®P'^®sents the step size of the Gj^ grid, We 

can represent the iteration operator (e.g., eq. (13))on the finest grid 

G„a5; 

L„ OCm) = fM (15) 

where i(i|n is the exact discrete solution on the G„ grid. We can write 

^ ''m 

where Uj^j is the approximate solution and is the error. Then we have the 
residual equation : 

(''«> = fM ■ l-M (“m> 

where R[^ is the residual of the approximation Uj^ on the Gj^ grid. "Lf.^ is in 
general different from in the nonlinear case, which complicates matters. 
Nevertheless, if is smooth, the error will be smooth, and the residual 
equation (17) can be solved on a coarser grid. Thus, for example, we can 
write 

f 

*^M-1 ^M ^ 

where w is an approximation to the error V|^ on the G|^_^ grid, and denotes 
interpolation from the G|^ to G^. After solving the problem (18) (usually with 
homogeneous boundary conditions), we interpolate the function W|^_i back onto 
the Gj^ mesh, and thus form an improved approximation: 

^“m^ new " old ^M-1 ^^M-1^ 


In the complete MG algorithm, the solution of equation (18) is also 


performed by relaxation; and if the convergence rate falls below a prescribed 
level, we can apply a similar procedure, backing up to the level, 

and sc on, until we arrive at G^, if necessary . The G^ grid is so coarse 
that a direct solution could be used economically, but we have used iteration 
here also. 

Full Approximation - In the general nonlinear case, the form of the 
operator L can be quite complicated — more so than the original operator, 

L — and thus applications to, say, the full potential equation may be tedious 
to program. It turns out that for the transonic small disturbance equation, 
the job is simple, and our first program did use the exact expression for L 
in an efficient way. However, there is an equivalent, easier method for 
solving the residual equation, which we call the full approximation method, 
as follows: • * 

Suppose we add to both sides of equation (18) the function 

•-M-l “ ^M-1 '^M-1 (21) 


Then, since (Wf^_^) + (u^^) 

we have 

'*M-1 " ^M-1 " ^ (22) 

We can now use the original operator on all the grids, which greatly simplifies 
the programming. The right-hand side of equation (22) is the difference 
between the residuals of u^^ calculated with the coarse-and fine-grid operators. 


Note that when the solution converges on the G|^ grid, then 

- 0 

(Rm) - 0 


0 


7 


(23a) 

(23b) 


1 


but win remain finite, since (jj^ is a solution on the G|^ grid; R,^_i is 
essentially the truncation error of the operator, 

After equation (22) is solved to sufficient accuracy, we determine the 
function 

Wm.! “ ‘J'm ” ^ 

by subtraction at all points of the grid and then interpolate to 
the grid as before in equation (19). 

More explicit details of the method will be deferred to a forthcoming 
report. 

RESULTS AND DISCUSSION 


In order to estimate the efficiency of the method, a work unit can 
be defined as the amount of computational effort ‘required for one relaxation 
sweep on the (finest) grid. Thus a relaxation sweep on the Gj^ grid costs 
n^y = (1/4) work units, for example, Likewise, when we calculate tl,e residuals 
for the G|^ grid, we perform these calculations at the points of the grid, 
i.e., 1/4 as few points; hence each residual calculation costs less than 1/4 
the effort of a relaxation sweep on the Gj^ grid, or approximately (1/4)^'*^’*'^. 

Note that this is an overestimate, since the tridiagonal system (13) is not 
inverted, nor do we calculate the iteration coefficients during the residual 
calculations* On the other hand we did not count the work of interpolation 
in equation (19), for example, or any other "overhead" of that type. 

An overall estimate of efficiency can be given by the number 



I 


.X— L 


where 


lj|“ norm of after first sweep on 
||R(^, I* norm of Rj^ after work units 


and 


Rj^ * AyS^jRj^) 


1/2 


( 26 ) 


Hence the norm we use is the root mean square of the residual on Gj^, This 
number is typically about 5 to 10 times smaller than ^he maximum norm in 
transonic problems. We consider an approximate solution to be converged 
when 

||R|^^|| < C/(no, of grid points) (27) 

where the prescribed constant C is typieally chosen as i so as to estimate 

the nominal truncation error. - ' 

Unstretched Grids 

In the case of a grid with constant steps in both directions, the present 
MG method performed quite well. In the following some typical results are 
summarized. All of the figures shown are copies of the screen display, on 
a remote computer terminal, of an abbreviated history of the MG runs. The 
first integer is the grid level, M, corresponding to Gj^ in our text. The 
next three "E"-format numbers are: 


1 . max 
ij 




2 . 


M 


3. max 
id 


Tij 


(See equation (13)). 
(See equation (26)). 
(See equation (14)) . 


The two integers following 1. above are the i,j location where the maxlR.. J 

id 
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occurred. The last two numbers 1n a row are the number work units, n^, 

W 


and the number of supersonic points. One row is printed for each relaxation 
sweep on the finest (Gj^) grid, but not for the coarser grids. However, each 
time the calculation "backs up" to a coarser grid, the words RESCAL are 


printed and the value of max 

id 


R.. and 


M 


are printed, 


together with the grid level (L) which has just been relaxed. Note that these 
norrrts correspond to in the right-hand s'ide of equation (22). In 

all MG runs shown, a relaxation factor of 1.0 was used on all grids. Like- 
wise all MG runs in these examples used five levels of grids (M « 5), with 
G^ being 4x2 mesh cells in the x-and y-directions, respectively, and 6g 
being a 64 x 32, We have done 6 1<2/els, with Gg being 128 x 64 with no 
deterioration in MG performance. 

Laplace*s Equation . - To show just how fast the MG method works for 
a nice, smooth, elliptic problem, we present a run in figures 1 and 2 for 
Laplace's equation with the prescribed normal derivative equal to sin ttX 
along y = 0. In figure 1, the convergence history is shown for G^, a 32 x 16 gri 
and according to equation (23), we achieved a = .540. Now because of the 
smoothness of the solution, it may be expected that interpolating the con- 
verged solution onto Gg will give a very good starting approximation for 
Gg. This is true, as can be seen on figure 2, where the Gg grid was started 
with the interpolated G^ solution. For Gg, we obtained a >= .583, but the 
efficiency of the two combined levels is more like a = .461 

Figure 3 shows the convergence history for the same problem, using 
SLOR all the way on Gg, achieving convergence in n^ = 141, yielding a = .924. 
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Nonlinear, Subcritical Flow « 0.7). - In figure 4 Is shown the 
history of a nonlinear, but subcritical flow solved by MG. Here the con- 
vergence rate Is the same or better, on 6g, as It was for Laplace’s equation 
with smooth boundary conditions discussed previously (l.e,, a * .549). In 
this case, however, the Neumann boundary condition Is an ”N-wave" - far from 
smooth - and hence we can conclude that discontinuous boundary conditions do not 
deteriorate MG performance. SLOR, with « » 1.85, achieved a » .868. 

Supercritical Flow (M„ « .85). - Figure 5 Illustrates the history for 

i I 

a typical supercritical flow with a moderate-sized supersonic region. Since 
the Gg grid has 2145 grid points, the supersonic region, with 124 points, 
occupies about 6 % of the grid. For this case, a - .593. The same case 
using SLOR all the way converged In « 58, using a relaxation factor of 
1.85, and the a = .855. , . 

Highly Supercritical Flow (M^ « ,95). - Figure 6 illustrates the 
history for a highly supercritical flow, where the shock wave is at the 
trailing edge of the airfoil. Note the final number of supersonic points 
(355) is established after 38 work units. It is typical that, at that point, 
the MG method begins to work best, since most of the high frequency error 
components have been eliminated, For this case, a - .858, achieving con- 
vergence in 67.6 work units. The same case was converged with SLOR all the 
way in 228 work units, with a = .957. 

Stretched Grids 

We found quickly that vertical line relaxation alone is not the best 
way to relax the solution in the MG mode in the case of a stretched grid. A 
possible explanation for this is that all of the high-frequency error 


components are not rapidly damped by vertical line relaxation in a general 
stretched mesh, where the mesh aspect ratio varies from very small to very 
large values. For if we consider the line relaxation algorithm for Laplace's 
equation, with a local mesh aspect ratio equal to A, and a relaxation factor 
w, the amplification factor isj 


A [2(l-w) +0)6^®^] 

^ ^ A(2-fae + 2w (1-cos 0y) 


(as) 


If A “(Ay/Ax)^ is large, we have a problem, for then, with » 0 and Oy » — , 


we have g(0, A,w) « 

c. A(2°m) 4 


+ 2m 


m 


and if w ® 1 , we see that 


g 


as A — « 


Clearly, choosing w « 2 alleviates the problem, but then other high-frequency 
components are retarded, i.e., for and Oy = 0, 


g( I , 0; A, w) 


/4(l-u)^ + J- 


4 + m 


2 


(30 


which approaches 1 as w nears 2.0. A solution to this problem is to sweep 
in all directions alternately (forward, backward, up, and down, in a general 
problem), but of course special care must be taken in supersonic regions. 

Figure 7 shows an M6 run with vertical line relaxation for the = .95 
flow, with the grid stretched to infinity in both the x- and y-directions. 

A logarithmic stretch was used, with 30, j of the grid points in the x-direction 
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on the airfoil chord. Note that the maximum residual tends to occur far 
above the airfoil (small values of j), where Ay/ Ax is large. For this case, 

a « .936. The same case, solved by SLOR takes about 382 cycles to 
converge (a = .974). Some benefit is still achieved from the MG 

mode of operation; even though the MG performance is far worse than what we 
believe can be obtained by a better relaxation algorithm. 

Since this last case is a particularly interesting flow, we have 
included some pictures of t%e output for the pressure distribution along 
y = 0 (Fig. 8), a chart of the Mach numbers in the computational plane 
(Fig. 9) and an isobar plot (Fig. 10). Note in figure 8 that &n obi ique 
shock occurs at the trailing edge, followed by a nearly-constant velocity 
supersonic zone in the wake, then a normal shock in the wake about 1 /2-chord 
behind the trailing edge, and finally a very, slow recovery to free-stream 
conditions. The airfoil lies between I = 24 and 42 ( x < .5). Figures 9 
and 10 show the "fishtail" shock pattern more clearly. In figure 9, only 
odd values of J are printed in order to fit the picture on the screen. 

J = 1 corresponds to infinity, as do I = 1 and 65; The values of I are the 
first column of integers, and the Mach numbers x 100 are shown in the array. 
Flow Is from top to bottom in the picture, with the line y = 0 (and the 
airfoil surface) on the left (J = 33, see bottom row of integers indicating 
the value of J). The isobar plot, figure 10, uses integers for supersonic 
flow values. The triangular region of nearly-constant velocity between the 
oblique shock at the trailing edge' and the normal shock in the wake is 
clearly evident. 

A summary of all these results is shown in figure 11. 
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CONCLUDING REIARKS 


The inultigrid (MG) method for accelerating relaxation calculations 
has been shown to be applicable to transonic flow with embedded shock waves. 

In this paper, vertical line relaxation was used for solving the nonlinear, 
conservative difference equation modelling the small-disturbance equation for 
the velocity potential. The multigrid approach appears to work about three 
to five times faster than optimal SLOR on unstretched grids of moderate 
size (64 X 32); The relative advantage of MG to SLOR increases as the grid 
gets finer, since the MG convergence rate is nearly independent of mesh size. 

On stretched grids, the present MG method slows down, being only about 
twice as fast as SLOR, It is felt that the reason for this is clear; 
the indicated remedy being alternating-direction relaxation sweeps. 

Future investigations will include the alternating sweeps, and the 
extension of the method to lifting flows. 
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APPENDIX 








Iteration Coefficients 

We have used various choices for iteration coefficients in equation (13). 
The coefficients used to make the calculations presented in this paper are 
simply based on the Newton linearization of equations (8), (9), and (11). 

They are as follows: 

First define: (dropping the j index, since all quantities are 

evaluated at the same i) 


I K- (y+DM^ (4>i+i - 
^ ^ Ax 



(Al) 



I 


I 


I 

I 

i 


\ 

i 


Then we have 

u, p = U,-4x-2 

A = C= -4y'^ 

B = + 2 0-p,) U^/u - n-., b,_ 1 

D = (1-M,) b,_ I - 2 m, . I U,., 

E = Mi-1 b,. 3 

Where = 0 if U^->0 

=1 u.<o 


(A2) 

(A3) 

(A4) 

(A5) 

(A6) 

(A7I 
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Figure 1. - MG solution of Laplace's equation with smooth boundary conditions. 32 x 16 grid. 
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Figure 3. - SLOR solution of Laplace's equation with smooth boundary conditions. 64 x 32 grid. 
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Figure 4. - MG solution of parabolic-arc airfoil , = 0.7, t = 0.1, 64 x 32 grid. 
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Figure 5. - MG solution of parabolic-arc airfoil. - 0.85, x - 0.1, 64 x 32 grid. 
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Figure 6. - MG solution of parabolic-arc airfoil. M = 0.95, t = 0.1, 64 x 32 grid 
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Figure 7. - NG solution of parabolic-arc airfoil. — 0.95, x = 0.1, 64 x 32 stretched grid 
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Figure 7. - Concluded 
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Figure 8. - Pressure distribution along y = 0 for solution of figure 7. 
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Figure 9. - Mach number chart for solution of figure 7 
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Figure 10. - Isobar chart for solution of figure 7 
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